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■ We investigate the notion of Bose-Einstein condensation of interacting particles. The 
' definition of the condensate is based on the existence of the dominant eigenvalue of the 

^ , single-particle density matrix. The statistical properties and the characteristic temperature 

^ ■ are computed exactly in the soluble models of two interacting atoms. 

Q' 

(N 

m '. I. INTRODUCTION 

> 

ly-^ ' The recent observation of the quantum degeneracy in mesoscopic samples of dilute alkaline, hydrogen and 

[ metastable helium gases marks an important breakthrough in a long standing quest for the experimental 

. realization of the Bose-Einstein condensation in its original form |^,||. Compared to other systems - cold 

' gases, say, or liquid helium - these new systems are characterized by an unprecedented range in life-time, 

, density, and temperature, by their purity and variability, and by the precision with which they can be ma- 

nipulated and controlled. The precision control of the system parameters, including the interaction strength 
Q and the particle number is particularly remarkable, as it offers for the first time the opportunity to 

e study in greater detail fundamental issues, like the transition from the non- interacting gas to the interacting 
I I gas, the emergence of a condensate in going from few to many particles, or the influence of the interaction 

. on the structure and the dynamics of the condensed phase. 

' It is easy to understand the notion of Bose-Einstein condensation of the ideal gas. In the non-interacting 

O ■ system the total energy is just the sum of energies of each particle, each particle having at its disposal the set 

O ■ of single-particle energy states. Compared with their distinguishable counterparts, indistinguishable Bose 

particles have a pronounced tendency to bunch in the same state Thus, as temperature is lowered, and 
the total energy becomes limited a resource, condensation occurs, which is the occupation of the single- 
. ^ particle ground state by a macroscopic fraction of all particles. Of course, for ultra low temperatures also 

■ distinguishable particles would pile up in that state. But for Bose particles, as the temperature is lowered, 
condensation occurs suddenly, at a well defined temperature, and none but only the single-particle ground 
state displays a truly macroscopic population in the low temperature regime. 

The picture is not so clear for interacting Bose particles. Collisions certainly contravene the Bose particles' 
tendency to pile up in a single state. In contrast to the ideal Bose gas, no exact solution is known for the 
iV-particle Hamiltonian describing the particles in the external trapping potential (to the delight of theorists 
it is a nearly perfect harmonic oscillator potential) even if the particles' pair interaction is modeled by 
a zero range pseudo-potential. Even simpler a problem, that of a gas in a three-dimensional box with 
periodic boundary conditions, is still poorly understood beyond the level of mean-field description. Many 
different predictions for the interaction induced shift of the condensation temperature, for example, coexist 
in literature Q. 

The minimal configuration where indistinguishability and pair-interaction become manifest is a "gas" of 
two particles. Hence, a two-particle system is the ideal starting point for a thorough investigation of the 
role of indistinguishability, the mechanism of the BEC, and its dependence on the pair interaction strength, 
without invoking the mean field theory approach. 

Our discussion will be based on two distinct models. The first model is due to Elliot H. Lieb It 
describes particles traveling on a circumference of a circle with repulsive (5-like interaction. The exact wave 
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functions for the N atom problem are analytically known in this case, while the eigenvalues of the energy are 
computed numerically from a finite set of transcendental equations. Of course this being a one-dimensional 
model, it does not display a phase transition in the thermodynamic limit But since we are going to 
study only a case of a finite number of atoms, no non-analyticity of any kind is expected anyway. In fact, 
for the purpose of illustration we restrict our attention to the system of two atoms only. This will prove to 
be sufficient to illustrate the general strategy, which should be used in experimentally more relevant cases. 
The direct extension of our calculations in the Lieb model to more than two atoms is very tedious but is in 
principle possible. The second model is that of two atoms in a one-dimensional harmonic oscillator trap. In 
this case the extension to three dimensions (3D) is simple but going beyond 2 atoms is very difficult. 
In this model, the exact eigenstates of the N = 2 particle Hamiltonian are expressible in terms of known 
special functions and the energy eigenvalues again require a numerical solution of a transcendental equation. 

Our paper is organized as follows: In section || we present a general approach to the Bose-Einstcin 
condensation of a finite interacting system, which is based on the properties of the single-particle density 
matrix. In section [II this general strategy is applied to the Lieb model, while in section IV we apply this 
strategy to the model of two atoms in a one-dimensional (ID) harmonic oscillator trap. The definition of a 
critical temperature of our two models is then examined in section The two final sections are devoted to 
concluding remarks and acknowledgements. 



II. INTERACTION 



Condensation in an ideal gas refers to the macroscopic occupation of the single-particle ground state. But 
what is the condensate for the interacting Bose gas? One may be tempted to call the multi-particle ground 
state of the system the condensate, but this is wrong. To repeat - the condensate is characterized by (i) a 
single-particle state which becomes (ii) macroscopically populated in the low temperature regime. Which 
state could that be? Of course, the correct definition should reduce to the one for the ideal Bose gas when 
the interaction strength tends to zero. The correct definition, it turns out, was already implicit in the study 
of the coherence properties of the superfluid helium, where the neutron scattering, probing the single-particle 
density matrix has revealed the presence of an off-diagonal long-range order |ll| . 

The single-particle density matrix derives from a generalization of the full density matrix, gi^^^ = Tr'[£)], 
where g is the full density matrix of the A^-particle system and Tr' refers to a trace over the degrees of 
freedom of a subset of — 1 particles. In the coordinate representation 

g''^\x,x')^ J dx2---J dx]yg{x,X2, ■ ■ ■ ,xn;x' ,X2, ■ ■ ■ ,xn) ■ (1) 

As a Hermitian, positive operator with a unit trace, the single-particle density matrix has a spectral decom- 
position 

g^'\x,x') = Y,X,cf>*{x)cf>,{x'), (2) 
j 

which defines a complete set of single-particle states and corresponding eigenvalues Aj > 0, J^j -^j — ^■ 

The j-th eigenvalue measures the weight of the corresponding single-particle state (often called a natural 
orbital) in the spectral decomposition. It may be represented as Xj — -^-j^, where (uj) may be identified with 
the mean number of particles in the j-th single-particle state. We are entitled to call one of the states, (pQ 
say, condensed if the corresponding eigenvalue is of order 0(1), all the other eigenvalues being much smaller, 
Aj < Aq. 

To compute the single-particle density matrix and its ensuing spectral decomposition is not an easy task 
in the general case, where the full density matrix is given in terms of the system's statistical specifications. 
For a iV-particle system in contact with a thermostat at temperature T, for example, the density matrix 
reads g = ^e~^^ , where /3 = l/(fcBT), H is the system Hamiltonian, and Z = Tig is the canonical partition 
function. 

Being a function of the system Hamiltionian, the density matrix is diagonal in the energy representation, 
g{xi,...,XN;x[,...,x'j^) = ^^e~''^^'■lp^{xl, . . . ,xn)iI'i{x[, . . . ,x'n) , (3) 

i 

where the energy values Ei and associated multi-particle wave functions ipi are given by the solution of the 
stationary Schrodinger equation 
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Hip{xi, . . . ,xn) = EiIj{xi, . . . ,xn) ■ (4) 

In the case of non-interacting particles the only correlation of the multi-particle density matrix comes 
from bosonic symmetrization. Hence the natural orbitals of the single-particle reduced density matrix (^) 
are simply the eigenstates of the single-particle Hamiltonian. 

In the case of interacting particles the situation is more complicated. Yet the symmetry may be of 
help. For the system in a box with the periodic boundary conditions, the equilibrium single-particle density 
matrix must be (i) invariant under spatial translations, and (ii) it must be spatially periodic. The spectral 
decomposition is therefore realized by the Fourier expansion and the natural orbitals are just plane waves 
or the momentum states, 

p( >{r,r )^}_^— (5) 

p 

where the momentum p — ^{ji, j2, js) (ji J2, and are integer numbers). The situation is very analogous 
in the standard theory of homogeneous Bose systems which explains why the condensate is identified with 
the zero momentum state. 

In the inhomogeneous case, relevant for the present experiments with magnetically trapped gases, symme- 
try alone cannot determine the natural orbitals. In the most regular case of spherically symmetric trap, the 
natural orbitals may be chosen proportional to the spherical harmonics (the density matrix may be diago- 
nalized together with the angular momentum operators), but their radial dependence feels the interaction 
and in general changes with temperature. 



III. THE LIEB MODEL 



A model which to a large extent allows for an explicit construction of the single-particle density matrix 
is the Lieb model The Lieb model describes a collection of interacting particles on a circle, the pair- 
interaction potential being modeled by a repulsive (5-function. Although this choice of an interaction potential 
may seem strange it does describe very well collisions between very slow particles, which in a gaseous phase 
dominate at low temperatures, and has been widely used in the context of BEC 

We commence with two free particles subject to the aforementioned conditions. We introduce conveniently 
scaled variables that will be used hereafter: the circle's circumference L defines the unit of length and 

/2mL^ defines the unit of energy. In scaled units the Schrodinger equation for two free particles reads 



dx\ dx2 



1p{xi,X2) = E'lp{xi,X2) , (6) 



where Tp{xi,X2) must obey Bose-symmetry 

■4>{X1,X2) = 1p{x2,Xl) , (7) 

and periodic boundary conditions, 

ip{xi + n,X2 + m) = ip{xi,X2) , rt, m integer . (8) 
Solutions of equations (0)-® are given in terms of plane waves 



which posses kinetic energy 



E = kl + kl. (10) 

Due to the boundary conditions ^ the wave numbers are quantized 

fc, = 27rz/,, I/, =0,±1,±2,... . (11) 

As for any pair of wave numbers {ki,k2) the wave functions ijjkik2 and "(Afegfei I'efer to the very same state, 
we impose the restriction 
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ki < k2 



(12) 



in order to avoid the double counting of states. 

We now turn to the problem of two interacting particles. The Schrodinger equation pertinent to this 
problem reads: 
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dxl 



+ 2c5{x2 - Xi) \ tl){xi,X2) = E'iIj{xi,X2) , 



(13) 



where c is the strength of the interaction. Safe for the collisions, which occur for Xi = X2, we just face the 
Schrodinger equation for two free particles. Across the line xi = X2 the first spatial derivatives of the wave 
function displays a step discontinuity (because the derivative of the step function is the 5 - function) , 



^ a_ 

dx2 dxi 



_d a_ 

8x2 dxi 



2CV' 1x2=2:1 • 



It is sufficient to find a solution in the region 

S : < xi < 2:2 < 1 



(14) 



(15) 



which is a simplex in the {xi, X2) plane, for if we know '4'{xi, X2) in S, we know it everywhere by symmetric, 
suitable for bosons, extension. 

In the simplex S the problem reduces to solving the free Schrodinger equation (^ supplemented by the 
boundary conditions 



(16) 



d d 

^{0,X2) = '(p{x2,l) -Q^ik{x,X2)\x=a = ^■4^{x2,x)\x=i . (17) 

Here equation (|l6| ) is the adaptation of ( |l4[ ) for the simplex S*, and equation (|l^) is the adaptation of the 
periodic boundary conditions (||). 

One should expect the solution to the problem of two interacting particles to be still a combination of plane 
waves (as in the free case). After all, the atoms move freely most of the time - the only perturbation taking 
place when they collide. The particular form of the solution, which indeed respects that expectation, was 
proposed in a classic paper by Elliott H. Lieb and W. Liniger in 1963 ||] (where the authors dealt with the 
problem of N interacting atoms; to date, this so called Lieb model is the only soluble model of N interacting 
bosons with short-range interactions). Recently, intensively studied has been a model of interacting bosons 
with the interatomic potential proportional to the product of two particles' coordinates In this case the 
total Hamiltonian is just a quadratic form. 

In a slight deviation from the Lieb ansatz, the solution is constructed with the help of a permanent of 
plane waves, 

r7(xi, X2) = Per |exp(iA:iXj)| . (18) 

which is just like a determinant, yet with all minus signs turned into plus signs. In terms of the permanent 
(pisl), the solution reads 

V'(xi,X2) oc ^1 + ^— — ^ ?7(a;i,a;2) . (19) 

with di = d/dxi. As the solution appears as a combination of plane waves, the energy eigenvalues assume 
the same form as in the free case, see (^0|). 

The solution is specified by a pair of quasi- wave numbers (fci, fc2), the possible values of which are de- 
termined by the boundary conditions. Inserting the ansatz ( |l9|) into (|l^ ) and ([l7|), these conditions turn 
into 

K = TT^, ^ = 0,±1,±2,... , (20) 

A: = 2arctan f-^) + TTi^ , v = 0,1,2, ... . (21) 
\2k/ 
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where K = {ki + k2)/2 and k ~ (^2 — ki)/2. Here, the first equation reflects the quantization of the total 
momentum, which is a conserved quantity. In the second equation we refer to the principal branch of the 
arcus-tangent, and we seek solutions with fc > to avoid double counting. 

There are two interesting limiting cases: weak interaction (c 0) and strong interaction (c — > oo). In 
the limit of weak interaction k « ^/c for v = 0, and k ttv + c/ {vk) for = 1, 2, . . .. These values connect 
smoothly with the free case k = ttv, v — 0, 1,2, . . .. Level repulsion is strongest for synchronized motion, 
i. e. when both particles travel with approximately the same speed, fci ~ fc2: and thus have ample time to 
"feel" the collision. For the generic iV-particle system in a box, say, the level repulsion leads to a depletion of 
the density of low lying states, and the concomitant appearance of collective excitations with a phonon-like 
dispersion. 

In the limit of strong interactions k k, (v -\- l)7r for v <^ and k vk for 3> c. Low lying states of 
the strongly interacting two-particle Bose system, for example the ground state oc sin(7r(a;2 — can be 
identified with the two-particle states of an ideal Fermi gas. This is a particular instance of the Bose-Fermi 
transformation, which maps the one-dimensional hard-sphere Bose gas onto an ideal Fermi gas (see and 
references therein). 

In the general case the solutions for the quasi- wave number difference k may be enumerated k^, v = 
0, 1,2, . . ., the v-th solution being strictly bounded nv < k^ < tt^v -\- 1). Expressed in terms of sums and 
differences of quasi- wave numbers, the spectrum of energy eigenvalues is given by (see equation (p^) 

i;^, = 2fc2 + 27rV' , J^ = 0,1,..., M = 0,±1,.... (22) 

The corresponding eigenfunction may be also labeled, 'ipfj.iyi^iT ^2), where integer /i specifies the center-of- 
mass motion, and non-negative integer v specifies the relative motion. 

We are now in a position to study the system statistical properties. Adapting the defining equation (^) 
to our present situation, and considering for the time being only states with a zero total momentum, /x = 0, 
the system canonical density matrix is given by 



g{xi,X2;x[,X2) = ^p^g^{xi,X2;x[,x2) , (23) 
1^=0 

where pi, is the Boltzmann weight, 

P. = , (24) 

and Qi, is the pure state density matrix, 

Qi^ = 1pOu{xi,X2)lpo^{x[,X2) ■ 

It is instructive to compute the reduced single-particle density matrix g^f}}, see equation (|l]) for the defi- 
nition. After some tedious yet straightforward calculation one obtains 

(1). {-Y sin k^ - k„[l - {-Y cos k^]z 

) = fe. + (-).sinfc. '""'^'^'^ 

ki^ -\- ( — )^ sin ki, 

with z = l^; — The thermal equihbrium reduced density matrix is given by 

00 

g'^^\x,x') = ^p„gl^\x,x') . 
According to the discussion in section || it may be diagonalized 

C30 

g^^\x,x')= X,^,{x)^*{x'), 

j = -oo 

where \j = with (uj) the mean occupation of the single-particle state <pj{x). Due to the translational 
invariance of the Lieb model, the single-particle states are plane waves, (j)j{x) — e^'^*-'^, j = 0, ±1, ±2, . . .. 
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The state with the largest eigenvalue in this decomposition will be referred to as the single-particle ground 
state and denoted by 0o- In the present model (jjoix) = 1 - most simple a wave function indeed. 

We shall be interested in the mean and root-mean square fluctuations of the single-particle ground-state 
occupation, 

e= v'K)-(no)2. (26) 
The mean and mean-square occupation of the single-particle ground state can be expressed as 

(no) = 2 • P2 + 1 • Pi + • Po , (27) 



(nl) = 4 • P2 + 1 • Pi + • Po , 



(28) 



where P„ denotes probability to find n atoms in the single-particle ground state. 

The probability to find two particles in the single-particle ground state can of course not be calculated 
from the single-particle density matrix, but must rather be inferred from the full density matrix, P2 = 
{(j)o(f)o\g\(t>o4'o) ■ In the coordinate representation. 

In order to compute Pi we recall (tiq) = 2(0o|^?*•^'' l^o) by definition, which in combination with ( p7|) yields 

Pi 2 / dxdx' (f)l[x)Q'^^\x,x')<l)o{x') - 2P2 . 



Finally, as the statistical events to find zero, one or two particles in the single-particle ground state are 
mutually exclusive and complete, 

Po = 1 - Pi - P2 . 

We note that in a restricted ensemble with zero total momentum the probability to find exactly one particle 
in the single-particle ground state is strictly zero. Pi = 0. For if one particle would be at rest, the other 
particle must be at rest too. Hence, in such a restricted ensemble, only P2 and Pq can be different from zero. 

Before we continue with the description of our results, we must supply a few words about technicalities 
of this work. All of them regard computation of the thermal density matrix (^3|) . In order to compute the 
partition function Z we sum over the first 100 states. The number of states used in the construction of the 
thermal density matrix is chosen so that the sum of all probabilities Pi, ( p4[ ) differs from 1 by less than 0.01. 
This usually renders a number between v^ax = 1 (for very low temperatures) and Vmax = 6 (for relatively 
high temperatures, T sa 1000). This choice constitutes an important aspect as the time of computation 
grows very quickly with the number of states involved. 



<nn> 




20 40 60 80 
T [arbitrary units] 

FIG. 1. Temperature dependence of the mean ground-state occupation in the Lieb model for the free case(c = 0, 



squares) as well intermediate (c = 30, circles) and infinite interaction 



00, triangles). 
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Now we can analyze the dependence of the occupations, fluctuations and probabihties on temperature and 
interaction strength (defined by the magnitude of c) . Figure |l] presents the dependence of mean occupation 
of the ground state ( p7| ) on temperature for various values of the interaction magnitude c. Obviously, in the 
free case (c = 0) at T = both atoms are in the 1-particle ground state (which can be traced back to its 
definition) . For c 7^ mean occupation of the ground state at T = is still substantial, but certainly smaller 
than 2 and assumes the smallest possible value for c — > 00. Interestingly, for higher temperatures (T w 50) 
(no) can be greater in interacting cases than in the free one. Analogous behavior can be observed in the 
dependence of fluctuations (|2^) on temperature. Here only for the free case they can take on the zero value, 
but again they can be bigger than in interacting cases for higher temperatures (see figure ^). 

1.0 r 

w 0.8 
g 

0.6 : 

^ 0.4 ; 
0.2 ; 
0.0; 

20 40 60 80 100 

T [arbitrary units] 

FIG. 2. Temperature dependence of the fluctuations of the ground-state occupation in the Lieb model for the free 
case (c — 0, squares) as well intermediate (c = 30, circles) and infinite interaction (c 00, triangles). 




IV. TWO ATOMS IN A HARMONIC TRAP 



Before we turn our attention to an example of two interacting bosons in a harmonic trap, let us comment 
briefly on the non-interacting case. The Hamiltonian describing two non-interacting particles in a one- 
dimensional harmonic trap takes a very simple form: 

H ~ fiuj{a\ai + a\a2) , 

where ai and 02 are lowering operators of the first and the second particle, respectively. The canonical- 
ensemble density matrix is 



1 

exp 



Z 



H 
kf 



where the partition function Z ensures the normalization of the density matrix. Until now we have not 
specified whether we deal with distinguishable or indistinguishable particles. This last property determines 
the space of available states. For distinguishable particles the trace defining the partition function extends 
over all configurations of particles independently occupying the levels: 



Zei= J2 9"^+"^ = (i-9r' 



where q = exp[— Tiw/fcT]. For identical bosons there are fewer states since the quantum identical particles 
cannot be labeled. So the sum in the trace is restricted: 



zb=J2 9"' E = [(i-9)(i-9')]- 



712—0 ni— 712 

In both models discussed here we compute the probability of finding zero, one, and two particles in the 
"condensed" state. This is easily done for classical and quantum particles in the present example. The only 
tool needed is the sum of the geometric series. One easily gets 

pel _ 2 

Pf = 2q{l - q) 

P2' - (1 - 1? 
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for distinguishable particles and 



P? = g(i - e) 

Pi = {l-q){l-q^) 

for bosons. Knowing the probability distributions we may easily compute all expectation values. The 
simplest is just the temperature dependence of mean number of particles in the ground state. In the two 
cases considered we get 

(no)^' = 2(1 - q) 
(no)^ = (2-g)(l-g2) 

One can easily check that for all finite temperatures {uq)^ > (^o)^'- In other words, as the temperature 
decreases, the population of the ground state grows faster for bosons due to the restricted space of states. 
Our reader should have no problem extending the results presented so far to the case of N non-interacting 
particles in a ID harmonic trap. For N going to infinity this is the origin of the Bose- Einstein phase transition. 

We investigate now the behaviour of two identical, interacting bosonic atoms in a harmonic trap. Again 
we approximate atom-atom interactions with a repulsive potential of a zero range. This system has analytic 
solutions which have been derived in |p^ , where the three-dimensional case is discussed in greater detail. 
We summarize here the main results for the one-dimensional system. 

The Hamiltonian can be decomposed into two parts: 

P^ P^cin ^" P^rel 5 

where Hem is the Hamiltonian of the center of mass, whose coordinate is defined here as X — yjl/2{xi + 
a;2), and H^ei is the Hamiltonian of the relative coordinate, defined here as a; = yjl/2{xi ~ X2)- This 
unconventional choice of the two coordinates differs from the usual definitions by a numeric factor and 
makes the effective masses M of the two degrees of freedom equal. The Hamiltonian of the center of mass 
X is equal to the Hamiltonian of the harmonic oscillator, 

2M ^ 2 

whereas the Hamiltonian for the relative coordinate is 

H=-^^1 + Imcj'x' + 2cS{x). 

Here the last term, proportional to Dirac's (5-function, describes the zero-range interaction. The parameter 
c determines the intensity of the zero-range interaction. For c = 0, we have two free atoms. To introduce 

dimensionless coordinates we scale both x and X by the oscillator length d = y/ Now the Hamiltonian 
of the relative coordinate takes the following form: 

Hrel = nuj{-^W^ + ^X^ + 2cS{x)) , 

where the interaction strength depends on the dimensionless parameter c. Since we consider repulsive 
interactions, c must be greater than 0. 

The eigenvalues and eigenfunctions of the center of mass are well known: 

En,c7n = hLj{n+^) 

with n = 0, 1, 2, . . ., and 

All these eigenfunctions are symmetric under the exchange of the two particles. 

For the relative coordinate we have two different sets of solutions. Half of them are odd eigenfunctions 
which vanish for x = 0. Therefore they are not influenced by the zero-range interaction and reduce to the 
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usual solutions of the harmonic oscillator. We do not consider them, because they are odd functions of 
the relative coordinate and therefore their combination with any eigenstate of the center of mass describes 
fermions instead of bosons (incidentally, we note that two trapped fermions do not feel the zero-range 
interaction, in accordance with the well known fact that cold fermions are insensitive to s-wave scattering). 
The remaining eigenfunctions are: 

B t2 / 1 

where 



(here ip is the logarithmic derivative of the gamma function F) and U is the confluent hypergeometric function 
(see page 504 of reference |l4|). These eigenfunctions are symmetric under the exchange of the two atoms, 
so we can use them in conjunction with the eigenstates of the center of mass to describe two bosonic atoms. 
The values of energy are: 

En,rel = h.Uj{Vn + n = 0, 2, 4, . . . , 

where the non-integer parameters v are solutions of the transcendental equation 

cT{^v) = -T{l/2-u). 

They also satisfy the condition n < Vn < n + I and 

lim i>n — n (no interaction) , (29) 
c— >o 

lim M„ = n + 1 (infinite interaction). (30) 

c — *oo 

The eigenstates of the system are therefore: 

We can now build the density matrix for the two trapped atoms. At T = the density matrix is simply 
given by the ground state of the system 

and the one-particle density matrix is 

/oo 
dyp{x, y; x' , y) = ^ Xi(j)i{x)(j)^{x') (31) 
i=0 

The integration and the diagonalization as indicated in (^l]) must be done numerically. The single-particle 
ground-state wave function (po is plotted in figure || for different values of c. These states differ from the 
ground state of the harmonic oscillator more and more as the interaction increases. 




-3 -2 -1 1 2 3 4 5 

X (arbitrary units) 

FIG. 3. Single-particle ground-state wave functions of 2 atoms in a ID harmonic trap at T 
line), c = 1 (dotted line), c = 2 (dashed line), and c = oo (dash-dotted line). 



for c = (solid 
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For T > we work again in the framework of the canonical ensemble. The density matrix is 



exp[ 



T-3/2 



E 



V2 



' '2' 2 \ '2' 2 



(32) 



In this case, technical difficulties come from the numerical evaluation of the U -functions with large values of v, 
especially for moderate and strong interactions. This problem restricts our investigation to low temperatures 
(between and 6 nanokelvin, for a typical trap frequency w = 27r • 100 Hz). The partition function is given 
by the sum over the first 90 states, and again (like for the Lieb model) the density matrix is constructed 
so that the sum of probabilities Pnu = exp(— ^"^^ j^*^^ ) /Z differs from 1 by less than 0.01. The prescription 
for the calculation of occupations, fluctuations and probabilities Pq, Pi, P2 is the same as in the previous 
section. 

Note that, unlike in the translationally invariant case of the Lieb model, the natural orbitals (eigenvectors 
of the single-particle density matrix) depend here both on the temperature and on the interaction strength. 
The temperature dependence turns out to be relatively weak. In particular, the ground state (f>o of the 
thermal reduced density matrix is very similar to that at T — 0. This is consistent with the fact that the 
Gross-Pitaevskii equation ||l^ (used to determine the wave function of the weakly-interacting condensate 
within the so called mean- field approximation), which is obtained for T = 0, is a very good tool also for the 
conditions realized in the experiments, where T > inevitably. 




0.2 0.3 0.4 0.5 

q 

FIG. 4. Mean occupation of the single-particle ground state of 2 atoms in a ID harmonic trap vs. q = exp{—huj / kT) 
for c — (solid line), c = 1 (dotted line), c — 2 (dashed line), and c — 00 (dash-dotted line). 



The mean occupation of the single-particle ground state is plotted in figure In analogy to the Lieb 
model, (no) at T = reduces from 2 (no interaction) to lower values as c increases. However, the curves for 
different values of c do not cross each other. 




0.0 0.1 0.2 q 0.3 0.4 0.5 
FIG. 5. Fluctuations of the ground-state occupation for 2 atoms in a ID harmonic trap vs. q — exp( — Rtj/fcT) for 
c = (solid line), c — 1 (dotted line), c = 2 (dashed line), and c = 00 (dash-dotted line). 
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The fluctuations of the ground-state occupation are depicted in figure H The curves for different values 
of c seem to cross beyond the range of temperatures we have investigated. It is remarkable that for infinite 
interaction the fluctuations are a decreasing function of temperature, contrarily to the case of low c. 



V. THE CRITICAL TEMPERATURE 



The critical temperature Tc defines a range of temperatures Q < T < Tc where a macroscopic fraction of 
the gas occupies the single-particle ground state. It is therefore reasonable to define the critical temperature 
in our two- atom systems by means of the crossing between Pq and Pi or P2, in the same spirit of the approach 
presented in for a weakly interacting Bose gas. 



A. Tc for the Lieb model 



As regards the probabilities of populating the ground state by a certain number of atoms (Po,-Pi and P2), 
our first observation is that for any interaction strength there is a crossing between P2 and Pq for a non-zero 
temperature (see figure ^- the sample case of c = 30). It defines some characteristic temperature Tc below 
which the ground state is populated by 2 rather than no particles. We believe it is an analogue of the critical 
temperature for higher-dimensional systems in the thermodynamic limit. This crossing corresponds to the 
maximum of fluctuations = 1) as well as to the point where the mean occupation of the ground state 
(no) crosses this maximum becoming smaller than fluctuations for temperatures higher than Tc- All these 
features are illustrated in figure 0. 



2.0 




0.0 ^ ' ^ ' ^ ' ^ ' ^ ' ^ 

200 400 600 800 1000 

T [arbitrary units] 

FIG. 6. Temperature dependence of probabilities P2 and Po as well as of the mean and fluctuations of the 
ground-state occupation for the case c = 30 in the Lieb model. 



The most interesting aspect of the dependence of the characteristic temperature Tc on the interaction 
strength is whether it increases or decreases as interaction magnitude grows. This relation is presented in 
figure for all positive c values. We observe that initially the characteristic temperature increases with the 
interaction strength, but having reached a maximum value for c « 100 it decreases thus rendering some 
characteristic interaction magnitude. Such a maximum occurs also in real condensates |l7|. 




10 100 1000 10000100000 
c [arbitrary units] 

FIG. 7. Dependence of the characteristic temperature Tc on the interaction strength c - 
turn in the Lieb model. 



case of zero total momen- 
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We should remind ourselves the special condition of a zero total momentum imposed on the system. It is 
automatically fulfilled in the canonical ensemble for N oo. However, as here TV = 2 we should also look 
into the case when the total momentum is unconstrained. There are several important new features in this 
case - Pi is no longer zero and it may also cross Pq and P2- It is tempting to investigate two characteristic 
temperatures now: one defined by the crossing of P2 with Pq and the other one when Pi crosses Pq. Their 
dependence on the interaction magnitude c is presented in figure ^j. Circles indicate crossings between P2 
and Po whereas squares - between Pi and Pq. Now, as we see, for small c Pi - Pq is met first when going 
from high to lower temperatures. At some c, P2 - Pq temperature is higher. From some c between 10 and 
20 up there is no Pi - Pq crossing - simply Pq is always greater than Pi. If we follow the higher of the 
two temperatures as c grows, we notice that this quantity almost always decreases (except for very large c 
values) . This result is opposite to the case of zero total momentum. 



80 




0.1 1 10 100 1000 10000100000 



c [arbitrary units] 

FIG. 8. Dependence of the characteristic temperature Tc on the interaction strength c - case of any total momentum 
in the Lieb model. Circles indicate crossings between P2 and Po whereas squares - between Pi and Pq. 



B. Tc for the trapped atoms 



The correspondence between the maximum of fluctuations and the crossing between P2 and Pq that we 
have noted in the Lieb model might exist also for two bosons in a harmonic trap. Unfortunately it occurs 
at temperatures that are beyond the range of our investigation. Only for c = 00 and c = can this 
correspondence be safely confirmed. However, a qualitative extrapolation of our data (see figure ^ as an 
example) supports this conclusion. 




0.0 0.1 0.2 q 0.3 0.4 0.5 

FIG. 9. Fluctuations of the ground-state occupation (solid line) and probabilities P2, Pi and Pq for 2 atoms in a 
ID harmonic trap (dotted, dashed, and dash-dotted lines, respectively) vs. q — exp{—huj/kT) for c = 1. 



For two bosons in a trap the probability Pi docs not vanish at P 7^ 0, so the probabilities Pi and Pq 
cross at some characteristic temperature Td, although only for moderate interactions (c < 2). This critical 
temperature decreases from about 12 nK (c = 0) to about 3.8 nK (c = 2) for the used parameters. The 
crossing between P2 and Po occurs for any value of c and the critical temperature Tc2 goes from about 8.1 nK 
for free bosons to T = for c ^ 00 (see figure ^0|). We note the different behaviour of critical temperature 
when compared to the Lieb model. 
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0.5 1.0 1.5 2.0 

c (arbitrary units) 



2.5 



FIG. 10. Critical values Qc — exp[—huj / kTc] vs. c for 2 atoms in a ID harmonic trap. 



VI. CONCLUSIONS 



In this paper we have presented a detailed analysis of the statistical properties of two interacting bosons 
in one dimension. The principal advantage of our work lies in the fact that the general strategy could be 
carried out exactly in this case leaving no doubts on the nature of the approximations that would have to 
be employed otherwise. Still, the calculations involved do not present substantial mathematical difficulties 
permitting a complete review of the methods and results. 

The situation in a real experiment is much more complex. In recent experiments 2000 to 5 • 10^ atoms have 
been condensed in a nearly perfectly harmonic trap. While the statistical properties of the finite ideal Bose 
gas are well understood within the framework of the restricted ensembles (canonical |l^ and microcanonical 
|19| ), the impact of interactions is far from clear. Most phenomena are sufficiently well described within the 
mean-field approximation, suitable for the condensate at zero temperature. Finite temperature phenomena 
are significantly harder to compute. We have already mentioned conflicting predictions of the shift of the 
critical temperature . Another example is - as yet - unexplained temperature dependence of the measured 
eigenfrequencies of the condensate . 

We believe that a clear presentation of the fundamental concepts, applied to soluble models may shed 
some light on the problem of interacting cold bosons at a real-world scale. 
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